#!/usr/bin/env ruby
# coding: utf-8

#
# topological-hex, a program to compute combinatorial hexahedral meshes.
#
# Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
#
# List of the contributors to the development, description and complete License:
# see LICENSE file.
#
# This program (topological-hex) is free software:
# you can redistribute it and/or modify it under the terms
# of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program (see COPYING file).  If not,
# see <http://www.gnu.org/licenses/>.
#

gem 'slop', '~> 4.6'

require 'set'
require 'matrix'
require 'tempfile'
require 'pp'
require 'json'
require 'slop'

def reorder_quad(quad)
  min_i = quad.each_index.min_by { |i| quad[i] }
  Array.new(4) { |i| quad[(min_i + i) % 4] }
end

def unorient_quad(quad)
  quad = reorder_quad(quad)
  if quad[1] > quad[3]
    quad[1], quad[3] = quad[3], quad[1]
  end

  quad
end

HEX_NEIGHBORS = [[1, 3, 4], [2, 0, 5], [3, 1, 6], [0, 2, 7],
                 [7, 5, 0], [4, 6, 1], [5, 7, 2], [6, 4, 3]]

HEX_FACET = [
  [0,1,2,3], [4,5,6,7], [0,1,5,4],
  [1,2,6,5], [3,2,6,7], [0,3,7,4]
]

def unorient_hex(hex)
  min_i = hex.each_index.min_by { |i| hex[i] }

  neighbors_i = HEX_NEIGHBORS[min_i].sort_by { |i| hex[i] }

  adj_01 = HEX_FACET.find { |f|
    [min_i, neighbors_i[0], neighbors_i[1]].all? { |j| f.include? j }
  }.find { |i| i != min_i && i != neighbors_i[0] && i != neighbors_i[1]  }

  adj_02 = HEX_FACET.find { |f|
    [min_i, neighbors_i[0], neighbors_i[2]].all? { |j| f.include? j }
  }.find { |i| i != min_i && i != neighbors_i[0] && i != neighbors_i[2]  }

  adj_12 = HEX_FACET.find { |f|
    [min_i, neighbors_i[1], neighbors_i[2]].all? { |j| f.include? j }
  }.find { |i| i != min_i && i != neighbors_i[1] && i != neighbors_i[2]  }

  opposite = (0...8).find { |i|
    HEX_FACET.none? { |f| f.include?(i) && f.include?(min_i) }
  }

  [hex[min_i], hex[neighbors_i[0]], hex[adj_01], hex[neighbors_i[1]],
   hex[neighbors_i[2]], hex[adj_02], hex[opposite], hex[adj_12]]
end

def hex_faces(hex)
  [hex.values_at(0, 1, 2, 3),
   hex.values_at(4, 7, 6, 5),
   hex.values_at(0, 4, 5, 1),
   hex.values_at(1, 5, 6, 2),
   hex.values_at(3, 2, 6, 7),
   hex.values_at(0, 3, 7, 4)].map do |quad|
    reorder_quad(quad)
  end
end

def hex_corners(hex)
  [hex.values_at(0, 1, 3, 4),
   hex.values_at(1, 2, 3, 6),
   hex.values_at(0, 1, 2, 5),
   hex.values_at(0, 2, 3, 7),
   hex.values_at(5, 4, 6, 7),
   hex.values_at(5, 4, 7, 0),
   hex.values_at(4, 7, 6, 3),
   hex.values_at(7, 6, 5, 2)]
end

def flip_hex(hex)
  hex.values_at(0, 3, 2, 1, 4, 7, 6, 5)
end

class Mesh
  COMPONENT_ID = {'x' => 0, 'y' => 1, 'z' => 2}

  def self.load(file)
    File.open(file) do |io|
      load_from_io(io)
    end
  end

  def self.load_from_io(io)
    elements = []
    vertices = []

    while line = io.gets
      if line =~ /^\s*Hexahedra\b/
        n = io.gets.to_i
        elements.concat Array.new(n) {
          io.gets.split(" ")[0, 8].map { |i| i.to_i - 1 }
        }
      elsif line =~ /^\s*Tetrahedra\b/
        n = io.gets.to_i
        elements.concat Array.new(n) {
          io.gets.split(" ")[0, 4].map { |i| i.to_i - 1 }
        }
      elsif line =~ /^\s*Vertices\b/
        n = io.gets.to_i
        vertices.concat Array.new(n) {
          io.gets.split(" ")[0, 3].map(&:to_f)
        }
      end
    end

    new(vertices, elements)
  end

  def initialize(vertices, elements)
    @vertices = vertices
    @elements = elements
  end

  attr_reader :vertices, :elements

  def orient3d(*ids)
    matrix = Matrix.rows ids.map { |i|
      Vector.elements [1] + vertices[i]
    }

    matrix.determinant
  end

  def valid_corners?(hex)
    return false if orient3d(hex[0], hex[1], hex[3], hex[4]) < 0
    return false if orient3d(hex[1], hex[2], hex[3], hex[6]) < 0
    return false if orient3d(hex[0], hex[1], hex[2], hex[5]) < 0
    return false if orient3d(hex[0], hex[2], hex[3], hex[7]) < 0

    return false if orient3d(hex[5], hex[4], hex[6], hex[1]) < 0
    return false if orient3d(hex[5], hex[4], hex[7], hex[0]) < 0
    return false if orient3d(hex[4], hex[7], hex[6], hex[3]) < 0
    return false if orient3d(hex[7], hex[6], hex[5], hex[2]) < 0

    true
  end

  def num_valid
    elements.count { |hex| valid_corners?(hex) }
  end

  def fix_orientation
    return self if @elements.empty?

    elements = @elements.dup

    if num_valid < elements.size / 2
      elements.map! { |hex| flip_hex(hex) }
    end

    ordered_elements = [elements.pop]

    faces = Set.new
    hex_faces(ordered_elements[0]).each do |face|
      faces << face
    end

    until elements.empty?
      i = elements.find_index do |elt|
        faces.any? { |face| face.all? { |v| elt.include? v } }
      end

      elt = elements[i]
      elements.delete_at i

      new_faces = hex_faces(elt)

      if new_faces.any? { |x| faces.include? x }
        ordered_elements << flip_hex(elt)
        new_faces = hex_faces(ordered_elements.last)
        if new_faces.any? { |x| faces.include? x }
          raise "Invalid combinatorial orientation"
        end
      else
        ordered_elements << elt
      end

      faces.merge new_faces
    end

    Mesh.new(vertices, ordered_elements)
  end

  def boundary
    boundary = Set.new

    elements.each do |hex|
      hex_faces(hex).each do |face|
        face = unorient_quad(face)
        if boundary.include? face
          boundary.delete face
        else
          boundary << face
        end
      end
    end

    boundary
  end

  def boundary_vertices
    boundary.flat_map(&:itself).to_set
  end

  def changed_submesh(reference)
    old_set = reference.elements.map { |hex| unorient_hex(hex) }.to_set
    changed_elements = elements.reject { |hex| old_set.include? unorient_hex(hex) }
    self.class.new(vertices, changed_elements)
  end

  def to_s
    string = ""
    string << "MeshVersionFormatted 2\n"
    string << "Dimension 3\n"
    string << "\n"
    string << "Vertices\n"
    string << vertices.size.to_s << "\n"
    vertices.each do |v|
      string << v.join(" ") + " 0\n"
    end
    string << "\n"
    string << "Hexahedra\n"
    string << elements.size.to_s << "\n"
    elements.each do |hex|
      string << hex.map { |i| i + 1 }.join(" ") + " 0\n"
    end
    string << "\n"
    string << "End\n"

    string
  end
end

opts = Slop.parse do |o|
  o.string "-i", "--input",  "Mesh that needs to be modified"
  o.string "-o", "--output",
           "Output path for the renumbered mesh of the input file",
           :required => true
  o.string "-r", "--reference", "Reference mesh, numbered like the input mesh",
           :required => true
  o.string "-c", "--corrected", "Renumbered version of the reference mesh",
           :required => true

  o.on "-h", "--help", "Print this message" do
    $stderr.puts o
    exit 1
  end
end

a = Mesh.load(opts[:reference])
b = Mesh.load(opts[:corrected])

map = a.elements.zip(b.elements).each_with_object({}) do |(h_a, h_b), m|
  h_a.zip(h_b).each do |v_a, v_b|
    m[v_a] = v_b
  end
end

a = opts[:input] ? Mesh.load(opts[:input]) : Mesh.load_from_io($stdin)

a.vertices.replace b.vertices

a.elements.map! do |hex|
  hex.map { |v| map[v] }
end

puts a.boundary_vertices.map { |i| i + 1 }.join(", ")
File.write(opts[:output], a.to_s)
